RiverDrainage.f90 Source File
Define river network scheme
Source Code
!! Define river network scheme !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.2 - 22nd April 2024 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 10/Feb/2023 | Original code | ! | 1.1 | 18/Apr/2024 | ChannelInitiation modified for following DischargeRouting update | ! | 1.2 | 22/Apr/2024 | Flow direction convention imported from Morphology | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description ! routines to define river network scheme for discharge routing ! The original algorithms were initially implemented in a ! standalone program used to prepare input file to FEST. ! Starting from 2023 FEST release, river network ! topology orgamìnization is computed within the FEST ! before starting the simulation from the flow direction ! and flow accumulation maps. ! MODULE RiverDrainage !Modules used: USE GridLib, ONLY: & !Imported type definitions: grid_real, & grid_integer, & !Imported routines: NewGrid, & GridDestroy USE GridOperations, ONLY: & !Imported routines: IsOutOfGrid, & GetIJ, & GetXY, & CellArea USE DataTypeSizes, ONLY: & !Imported definitions: short, & float, & long USE LogLib, ONLY : & ! Imported Routines: Catch USE StringManipulation, ONLY: & !Imported routines: ToString USE Morphology, ONLY: & !Imported routines: HortonOrders , & CheckOutlet, & DownstreamCell, & DeriveSlope, & CellIsSpring USE GeoLib, ONLY: & !Imported definitions: Coordinate, & !Imported routines: Distance , & !Imported assignment: ASSIGNMENT (=) USE Utilities, ONLY: & !imported routines: GetUnit USE Morphology, ONLY : & !Imported variables: NW, & !North-West N, & !North NE, & !North-East W, & !West E, & !East SW, & !South-West S, & !South SE !South-East IMPLICIT NONE ! Local declarations: ! Local routines: PRIVATE :: SplitNetwork PRIVATE :: ExportReaches PRIVATE :: ExportShape TYPE Reach INTEGER :: id !==================================================== REAL(KIND = float) :: x0 !! beginning of reach x coordinate REAL(KIND = float) :: y0 !! beginning of reach y coordinate REAL(KIND = float) :: x1 !! end of reach x coordinate REAL(KIND = float) :: y1 !! end of reach y coordinate !==================================================== INTEGER :: i0 !! beginning of reach row INTEGER :: j0 !! beginning of reach column INTEGER :: i1 !! end of reach row INTEGER :: j1 !! end of reach column !==================================================== INTEGER :: ncells !!number of cells in a reach INTEGER :: order !! Horton-Sthraler order REAL(KIND = float) :: slope !! average reach slope (m/m) REAL(KIND = float) :: length !!reach length (m) REAL(KIND = float) :: area !! area drained by end section (m2) END TYPE Reach TYPE ReachNetwork TYPE(Reach),POINTER :: branch (:) INTEGER :: nreach END TYPE ReachNetwork TYPE ReachesList INTEGER :: id INTEGER :: i0 INTEGER :: j0 INTEGER :: i1 INTEGER :: j1 REAL (KIND = FLOAT) :: y0 REAL (KIND = FLOAT) :: x0 REAL (KIND = FLOAT) :: y1 REAL (KIND = FLOAT) :: x1 INTEGER :: n_cells REAL (KIND = FLOAT) :: slope ! (L/L) REAL (KIND = FLOAT) :: length !(m) REAL (KIND = FLOAT) :: area ! (cells) INTEGER :: order !Horton-Strhaler order TYPE (ReachesList), POINTER :: next !dynamic list END TYPE ReachesList !Global routines: PUBLIC :: BuildReachNetwork PUBLIC :: ChannelInitiation !private TYPE (ReachesList), POINTER, PRIVATE :: list, current, previous !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== ! Description: ! Define cells where channels begin. The routine searches for channel cells ! that have no other input channel cells. !SUBROUTINE MarkChannelBeginning ! !local declaration: !INTEGER :: i,j !INTEGER :: row, col !LOGICAL :: inputFound ! !-----------------------end of declaration------------------------------------- ! !DO i = 1, channelBeginning % idim ! DO j = 1, channelBeginning % jdim ! IF (channel % mat (i,j) /= channel % nodata) THEN ! inputFound = .FALSE. ! ! !check EAST cell ! row = i ! col = j + 1 ! IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN ! IF (flowDirection % mat (row,col) == W .AND. & ! channel % mat (row,col) /= channel % nodata ) THEN ! inputFound = .TRUE. ! channelBeginning % mat (i,j) = channelBeginning % nodata ! CYCLE ! END IF ! END IF ! ! !check SOUTH-EAST cell ! row = i + 1 ! col = j + 1 ! IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN ! IF (flowDirection % mat (row,col) == NW .AND. & ! channel % mat (row,col) /= channel % nodata ) THEN ! inputFound = .TRUE. ! channelBeginning % mat (i,j) = channelBeginning % nodata ! CYCLE ! END IF ! END IF ! ! !check SOUTH cell ! row = i + 1 ! col = j ! IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN ! IF (flowDirection % mat (row,col) == N .AND. & ! channel % mat (row,col) /= channel % nodata ) THEN ! inputFound = .TRUE. ! channelBeginning % mat (i,j) = channelBeginning % nodata ! CYCLE ! END IF ! END IF ! ! !check SOUTH-WEST cell ! row = i + 1 ! col = j - 1 ! IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN ! IF (flowDirection % mat (row,col) == NE .AND. & ! channel % mat (row,col) /= channel % nodata ) THEN ! inputFound = .TRUE. ! channelBeginning % mat (i,j) = channelBeginning % nodata ! CYCLE ! END IF ! END IF ! ! !check WEST cell ! row = i ! col = j - 1 ! IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN ! IF (flowDirection % mat (row,col) == E .AND. & ! channel % mat (row,col) /= channel % nodata ) THEN ! inputFound = .TRUE. ! channelBeginning % mat (i,j) = channelBeginning % nodata ! CYCLE ! END IF ! END IF ! ! !check NORTH-EAST cell ! row = i - 1 ! col = j - 1 ! IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN ! IF (flowDirection % mat (row,col) == SE .AND. & ! channel % mat (row,col) /= channel % nodata ) THEN ! inputFound = .TRUE. ! channelBeginning % mat (i,j) = channelBeginning % nodata ! CYCLE ! END IF ! END IF ! ! !check NORTH cell ! row = i - 1 ! col = j ! IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN ! IF (flowDirection % mat (row,col) == S .AND. & ! channel % mat (row,col) /= channel % nodata ) THEN ! inputFound = .TRUE. ! channelBeginning % mat (i,j) = channelBeginning % nodata ! CYCLE ! END IF ! END IF ! ! !check NORTH-EAST cell ! row = i - 1 ! col = j + 1 ! IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN ! IF (flowDirection % mat (row,col) == SW .AND. & ! channel % mat (row,col) /= channel % nodata ) THEN ! inputFound = .TRUE. ! channelBeginning % mat (i,j) = channelBeginning % nodata ! CYCLE ! END IF ! END IF ! ! IF ( .NOT. inputFound ) THEN ! channelBeginning % mat (i,j) = 1 ! END IF ! ! ELSE ! channelBeginning % mat (i,j) = channelBeginning % nodata ! END IF ! ! ! END DO !END DO ! ! !END SUBROUTINE MarkChannelBeginning !============================================================================== !| Description: ! Split channel network where a confluence of two different horton order ! channels occurs SUBROUTINE SplitNetwork & ! (orders, flowDirection, split) IMPLICIT NONE !Arguments with intent in: TYPE(grid_integer), INTENT(in) :: orders TYPE(grid_integer), INTENT(in) :: flowDirection !Arguments with intent out TYPE(grid_integer),INTENT(inout):: split !local declaration: INTEGER :: i,j INTEGER :: row, col INTEGER :: hortonOrder LOGICAL :: splitFound !-----------------------end of declaration------------------------------------- DO i = 1, split % idim DO j = 1, split % jdim !IF (channel % mat (i,j) /= channel % nodata) THEN hortonOrder = orders % mat (i,j) splitFound = .FALSE. !check EAST cell row = i col = j + 1 IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN IF (flowDirection % mat (row,col) == W .AND. & !channel % mat (row,col) /= channel % nodata .AND. & orders % mat (row,col) < hortonOrder) THEN splitFound = .TRUE. split % mat (i,j) = 1 CYCLE END IF END IF !check SOUTH-EAST cell row = i + 1 col = j + 1 IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN IF (flowDirection % mat (row,col) == NW .AND. & ! channel % mat (row,col) /= channel % nodata .AND. & orders % mat (row,col) < hortonOrder ) THEN splitFound = .TRUE. split % mat (i,j) = 1 CYCLE END IF END IF !check SOUTH cell row = i + 1 col = j IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN IF (flowDirection % mat (row,col) == N .AND. & ! channel % mat (row,col) /= channel % nodata .AND. & orders % mat (row,col) < hortonOrder ) THEN splitFound = .TRUE. split % mat (i,j) = 1 CYCLE END IF END IF !check SOUTH-WEST cell row = i + 1 col = j - 1 IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN IF (flowDirection % mat (row,col) == NE .AND. & ! channel % mat (row,col) /= channel % nodata .AND. & orders % mat (row,col) < hortonOrder ) THEN splitFound = .TRUE. split % mat (i,j) = 1 CYCLE END IF END IF !check WEST cell row = i col = j - 1 IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN IF (flowDirection % mat (row,col) == E .AND. & ! channel % mat (row,col) /= channel % nodata .AND. & orders % mat (row,col) < hortonOrder) THEN splitFound = .TRUE. split % mat (i,j) = 1 CYCLE END IF END IF !check NORTH-EAST cell row = i - 1 col = j - 1 IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN IF (flowDirection % mat (row,col) == SE .AND. & ! channel % mat (row,col) /= channel % nodata .AND. & orders % mat (row,col) < hortonOrder) THEN splitFound = .TRUE. split % mat (i,j) = 1 CYCLE END IF END IF !check NORTH cell row = i - 1 col = j IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN IF (flowDirection % mat (row,col) == S .AND. & ! channel % mat (row,col) /= channel % nodata .AND. & orders % mat (row,col) < hortonOrder ) THEN splitFound = .TRUE. split % mat (i,j) = 1 CYCLE END IF END IF !check NORTH-EAST cell row = i - 1 col = j + 1 IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN IF (flowDirection % mat (row,col) == SW .AND. & ! channel % mat (row,col) /= channel % nodata .AND. & orders % mat (row,col) < hortonOrder) THEN splitFound = .TRUE. split % mat (i,j) = 1 CYCLE END IF END IF IF ( .NOT. splitFound ) THEN split % mat (i,j) = split % nodata END IF !ELSE ! split % mat (i,j) = split % nodata !END IF END DO END DO END SUBROUTINE SplitNetwork !============================================================================== !| Description: ! Define channel cells. Two methods are possible: ! ! 1. area: channel begins where drainage basin exceedes a given area ! 2. ask: channel begins where AS<sup>k</sup> exceedes a given value, where A is ! the basin area (m<sup>2</sup>), S is the local slope (-), and k is the ! basin fractal dimension (default = 1.7) SUBROUTINE ChannelInitiation & ! (method, threshold, mask, flowAcc, flowDir, dem, channel ) IMPLICIT NONE !arguments with intent in: CHARACTER (LEN = *), INTENT(IN) :: method REAL (KIND = float), INTENT(IN) :: threshold TYPE (grid_integer), INTENT(IN) :: mask TYPE (grid_integer), INTENT(IN) :: flowAcc TYPE (grid_integer), INTENT(IN) :: flowDir TYPE (grid_real), INTENT(IN) :: dem !arguments with intent(inout): TYPE (grid_integer), INTENT(INOUT) :: channel !local declarations: INTEGER :: i,j,iu,ju,id,jd REAL :: area !! ( m<sup>2</sup> ) REAL :: scale = 1.7 REAL :: ask !!local value of AS<sup>k</sup> TYPE (grid_real) :: slope !----------------------end of declarations------------------------------------- !initialize channel to zero DO i = 1, dem % idim DO j = 1, dem % jdim IF (mask % mat (i,j) /= mask % nodata ) THEN channel % mat (i,j) = 0 END IF END DO END DO SELECT CASE (method) CASE ('area') DO i = 1, mask % idim DO j = 1, mask % jdim IF (mask % mat (i,j) /= mask % nodata ) THEN area = flowAcc % mat (i,j) * CellArea (flowAcc, i, j) IF (area >= threshold) THEN channel % mat (i,j) = 1 ELSE channel % mat (i,j) = 0 END IF END IF END DO END DO CASE ('ask') !compute slope in radians CALL DeriveSlope (dem, slope) DO i = 1, mask % idim DO j = 1, mask % jdim IF (mask % mat (i,j) /= mask % nodata ) THEN area = flowAcc % mat (i,j) * CellArea (flowAcc, i, j) ask = area * TAN (slope % mat (i,j)) ** scale IF (ask >= threshold) THEN !channel initiation found !set channel on all downstream cell iu = i ju = j DO channel % mat (iu,ju) = 1 CALL DownstreamCell (iu,ju,flowDir % mat(iu,ju),id,jd) IF ( CheckOutlet (iu,ju,id,jd,flowDir) ) EXIT iu = id ju = jd END DO END IF END IF END DO END DO CALL GridDestroy (slope) !other methods to be implemented END SELECT RETURN END SUBROUTINE ChannelInitiation !============================================================================== !| Description: ! Build reaches that compose the river network SUBROUTINE BuildReachNetwork & ! (maxReachLength, slopeCorrection, flowDirection, flowAcc, dem, & fileExport, shpExport, streamNetwork ) IMPLICIT NONE !arguments with intent(in): REAL (KIND = float), INTENT (IN) :: maxReachLength !!max length of a reach (m) REAL (KIND = float), INTENT (IN) :: slopeCorrection !! slope value to correct negative values TYPE (grid_integer), INTENT (IN) :: flowDirection !! flow direction TYPE (grid_integer), INTENT (IN) :: flowAcc !! flow accumulation TYPE (grid_real), INTENT (IN) :: dem !! digital elevation model INTEGER (KIND = short), INTENT (IN) :: fileExport, shpExport !arguments with intent (out): TYPE (ReachNetwork), INTENT (OUT) :: streamNetwork !local declarations LOGICAL :: outletSection LOGICAL :: endOfReach INTEGER :: countReaches INTEGER :: maxOrder INTEGER :: i,j,l !loop index INTEGER :: row,col !current cell INTEGER :: iDown, jDown !downstream cell INTEGER :: n_cells REAL :: reachLength TYPE (grid_integer) :: orderBeginning TYPE (grid_integer) :: orders TYPE (grid_integer) :: split TYPE (Coordinate) :: point1, point2 !- -----------------------------end of declarations---------------------------- !------------------------------------------------------------------------------ !(1.0) Initialize orders, orderBeginning and split grids !------------------------------------------------------------------------------ CALL NewGrid (orders, flowDirection) CALL NewGrid (orderBeginning, flowDirection) CALL NewGrid (split, flowDirection) !------------------------------------------------------------------------------ !(2.0) Build Horton orders grid !------------------------------------------------------------------------------ CALL HortonOrders (flowDirection, orders, maxOrder) !------------------------------------------------------------------------------ !(3.0) Build orderBeginning and split network !------------------------------------------------------------------------------ CALL MarkOrderBeginning (orders, flowDirection, orderBeginning) !split network CALL SplitNetwork (orders, flowDirection, split) !------------------------------------------------------------------------------ !(4.0) Build reaches !------------------------------------------------------------------------------ countReaches = 0 NULLIFY(list) ! at the beginning list is empty !initialize points to compute reach length point1 % system = flowDirection % grid_mapping point2 % system = flowDirection % grid_mapping DO l =1, maxOrder CALL Catch ('info', 'RiverDrainage', 'Elaborating reaches of stream order: ', & argument = ToString(l)) DO j=1,orderBeginning%jdim ! scan orderBeginning grid DO i=1,orderBeginning%idim IF(orderBeginning%mat(i,j) == l) THEN countReaches = countReaches + 1 reachLength = 0. n_cells = 0 row = i col = j endOfReach = .FALSE. outletSection = .FALSE. IF(.NOT.ASSOCIATED(list)) THEN ALLOCATE(list) !set first reach current => list ELSE ALLOCATE(current%next) !allocate next reach NULLIFY(current % next % next) current => current%next ENDIF current % i0 = i current % j0 = j current % order = l current % id = countReaches !follow the reach until an order change or to the outlet, !if it is the maximum order DO WHILE (.NOT. endOfReach) CALL DownstreamCell (row,col,flowDirection%mat(row,col), & iDown,jDown) n_cells = n_cells + 1 CALL GetXY (row, col, flowDirection, point1 % easting, point1 % northing) CALL GetXY (iDown, jDown, flowDirection, point2 % easting, point2 % northing) reachLength = reachLength + Distance (point1, point2) !length (m) IF ( CheckOutlet (row,col,iDown,jDown,flowDirection) ) THEN outletSection = .TRUE. endOfReach = .TRUE. ENDIF IF ( outletSection ) THEN current % i1 = row current % j1 = col current % n_cells = n_cells ! unit: cells (integer) current % length = reachLength ! !length (m) current % area = flowAcc % mat(row,col) ! * & !flowAccumulation % cellsize**2 ! if the reach length exceedes the maximum length => ! break the reach ELSEIF (maxReachLength > 0. .AND. & !if maximum length = 0 I do not break the reach reachLength >= maxReachLength ) THEN current % i1 = row current % j1 = col current % n_cells = n_cells current % length = reachLength !current % routingMethod = default current % area = flowAcc % mat(row,col) ! * & !flowAccumulation % cellsize**2 !(m2) !set elements for the following reach ALLOCATE(current % next) NULLIFY(current % next % next) current => current % next current % i0 = row current % j0 = col current % order = l countReaches = countReaches + 1 reachLength = 0. n_cells = 0 current % id = countReaches END IF IF ( .NOT. IsOutOfGrid(iDown,jDown,orders)) THEN IF (orders%mat(iDown,jDown) > l) THEN endOfReach = .TRUE. current % i1 = iDown current % j1 = jDown IF (n_cells == 0) THEN !this case occurs when a reach has just been terminated because too long current % n_cells = 1 CALL GetXY (row, col, flowDirection, point1 % easting, point1 % northing) CALL GetXY (iDown, jDown, flowDirection, point2 % easting, point2 % northing) current % length = Distance (point1, point2) ELSE current % n_cells = n_cells current % length = reachLength END IF current % area = flowAcc % mat(row,col) ! * & !flowAccumulation % cellsize**2 END IF END IF row = iDown col = jDown END DO ENDIF ENDDO ENDDO ENDDO !set other parameters !scan all list elements current => list DO !calculate slope from digital elevation model and check negative slope current % slope = ( dem % mat(current%i0,current%j0) - & dem % mat(current%i1,current%j1) ) / & current % length IF ( current % slope <= 0. ) THEN IF ( slopeCorrection > 0. ) current % slope = slopeCorrection ENDIF !calculate spatial coordinate CALL GetXY (current % i0, current % j0, dem, current % x0, current % y0) CALL GetXY (current % i1, current % j1, dem, current % x1, current % y1) IF ( .NOT. ASSOCIATED (current % next) )THEN !found last element of list EXIT END IF current => current % next END DO !------------------------------------------------------------------------------ !(5.0) remove temporary variables !------------------------------------------------------------------------------ CALL GridDestroy ( orders ) CALL GridDestroy ( orderBeginning ) !------------------------------------------------------------------------------ !(6.0) export files !------------------------------------------------------------------------------ IF ( fileExport > 0 ) THEN CALL ExportReaches ( ) END IF IF ( shpExport > 0 ) THEN CALL ExportShape ( dem, flowDirection ) END IF !------------------------------------------------------------------------------ !(7.0) populate stream network !------------------------------------------------------------------------------ ALLOCATE ( streamNetwork % branch (countReaches) ) streamNetwork % nreach = countReaches current => list DO i = 1, countReaches streamNetwork % branch (i) % id = current % id streamNetwork % branch (i) % i0 = current % i0 streamNetwork % branch (i) % j0 = current % j0 streamNetwork % branch (i) % i1 = current % i1 streamNetwork % branch (i) % j1 = current % j1 streamNetwork % branch (i) % x0 = current % x0 streamNetwork % branch (i) % y0 = current % y0 streamNetwork % branch (i) % x1 = current % x1 streamNetwork % branch (i) % y1 = current % y1 streamNetwork % branch (i) % ncells = current % n_cells streamNetwork % branch (i) % slope = current % slope streamNetwork % branch (i) % length = current % length streamNetwork % branch (i) % area = current % area streamNetwork % branch (i) % order = current % order current => current % next END DO RETURN END SUBROUTINE BuildReachNetwork !============================================================================== !| Description: ! Export reaches on file SUBROUTINE ExportReaches & ! ( ) IMPLICIT NONE !local declaration: INTEGER ( KIND = short ) :: fileunit !-----------------------end of declaration------------------------------------- fileunit = GetUnit () OPEN (unit = fileunit, file = 'stream_network.txt' ) WRITE (fileunit, '(a)') 'surface_routing_reaches' WRITE (fileunit, '(a)') ' ID X0 Y0 & X1 Y1 CELLS STHRALER_ORDER & SLOPE_L/L LENGTH_m DRAINED_CELLS ' !scan all list elements current => list DO WRITE(fileunit, '(I,3X,E,E,E,E,I,I,10x,E,1X,E,1X,E,I,11X,E,1X,E,& 1X,E,1X,E,1X,E,5X,E,E,E)') & current % id, current % x0, current % y0, & current % x1, current % y1, current % n_cells, & current % order, current % slope, current % length, & current % area IF ( .NOT. ASSOCIATED (current % next) ) THEN !found last element of list EXIT END IF current => current % next END DO CLOSE (fileunit) END SUBROUTINE ExportReaches !============================================================================== !| Description: ! Define cells where the beginning of a reach of different order takes place. SUBROUTINE MarkOrderBeginning & ! (orders, flowDirection, orderBeginning) IMPLICIT NONE !Arguments with intent in: TYPE(grid_integer),INTENT(in):: orders TYPE(grid_integer),INTENT(in):: flowDirection !Arguments with intent out TYPE(grid_integer),INTENT(inout):: orderBeginning LOGICAL :: outlet INTEGER :: row,col !current cell INTEGER :: iDown, jDown !downstream cell INTEGER :: reachOrder INTEGER :: maxPathOrder INTEGER :: i,j !- End of header -------------------------------------------------------------- DO j=1,orderBeginning%jdim DO i=1,orderBeginning%idim IF(CellIsSpring(i,j,flowDirection)) THEN !found a spring orderBeginning%mat(i,j) = 1 row = i col = j outlet = .FALSE. reachOrder = 1 DO WHILE (.NOT. outlet) ! follow the reach till the outlet CALL DownstreamCell(row,col,flowDirection%mat(row,col),iDown,jDown) outlet = CheckOutlet (row,col,iDown,jDown,flowDirection) IF (.NOT. outlet) THEN IF(orders%mat(iDown,jDown) == reachOrder + 1) THEN !found the beginning of a reach of increased order reachOrder = reachOrder + 1 orderBeginning%mat(iDown,jDown) = reachOrder ENDIF ENDIF row = iDown col = jDown END DO ENDIF ENDDO ENDDO !-----------remove duplicates------------ DO j=1,orderBeginning%jdim ! scan entire grid DO i=1,orderBeginning%idim IF(CellIsSpring(i,j,flowDirection)) THEN !found a spring row = i col = j outlet = .FALSE. maxPathOrder = 1 DO WHILE (.NOT. outlet) ! follow the reach till the outlet CALL DownstreamCell(row,col,flowDirection%mat(row,col),iDown,jDown) outlet = CheckOutlet (row,col,iDown,jDown,flowDirection) IF (.NOT. outlet) THEN IF(orderBeginning%mat(iDown,jDown) <= maxPathOrder) THEN !only one point where a N order reach begins is possible orderBeginning%mat(iDown,jDown) = 0 ELSEIF (orderBeginning%mat(iDown,jDown) > maxPathOrder) THEN maxPathOrder = orderBeginning%mat(iDown,jDown) ENDIF ENDIF row = iDown col = jDown END DO ENDIF ENDDO ENDDO RETURN END SUBROUTINE MarkOrderBeginning !============================================================================== !| Description: ! export shape file of river network SUBROUTINE ExportShape & ! ( dem, flowDirection ) IMPLICIT NONE !arguments with intent (in): TYPE (grid_real), INTENT (IN) :: dem TYPE (grid_integer), INTENT (IN) :: flowDirection !local declarations INTEGER :: K, jin,iin,ivalle,jvalle INTEGER :: dot REAL :: X,Y LOGICAL(4) :: res CHARACTER(300) :: command CHARACTER(300) :: filename INTEGER :: system !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ !(1.1) create "generate file" !------------------------------------------------------------------------------ OPEN(unit=876,file='tratti.gen') OPEN(unit=877,file='tratti.csv') WRITE(877,'(a)')'#ID;X0;Y0;X1;Y1;CELLS;STHRALER;SLOPE_L/L;LENGTH_m;& DRAINED_CELLS' !scan all list elements current => list DO WRITE(876,*) current % id jin = current % j0 iin = current % i0 WRITE(877,'(I8,";",E14.7,";",E14.7,";",E14.7,";",E14.7,";",I8,";",I8,";",E14.7,";",& E14.7,";",E14.7,";",I8,";",E14.7,";",E14.7,";",E14.7,";",E14.7,";",E14.7)') & current % id, current % x0, current % y0, & current % x1, current % y1, current % n_cells, & current % order, current % slope, current % length, & current % area DO WHILE (.not.((jin.EQ.current%j1) .AND. & (iin.EQ.current%i1))) X = DEM%xllcorner + jin * DEM%cellsize - DEM%cellsize / 2. Y = DEM%yllcorner + (DEM%idim - (iin-1)) * DEM%cellsize - DEM%cellsize / 2. WRITE(876,'(f20.7,a1,f20.7)') X,',',Y CALL DownstreamCell(iin,jin,flowDirection%mat(iin,jin),ivalle,jvalle) iin = ivalle jin = jvalle END DO !end single reach X = DEM%xllcorner + jin * DEM%cellsize - DEM%cellsize / 2. Y = DEM%yllcorner + (DEM%idim - (iin-1)) * DEM%cellsize - DEM%cellsize / 2. WRITE(876,'(f20.7,a1,f20.7)') X,',',Y WRITE(876,'(a3)') 'end' IF ( .NOT. ASSOCIATED (current % next) )THEN !found last element of list EXIT END IF current => current % next END DO !end cycle on all reaches WRITE(876,'(a3)') 'end' CLOSE(876) CLOSE(877) !!------------------------------------------------------------------------------ !!(1.2) create shape file !!------------------------------------------------------------------------------ filename = 'stream_network' !dot = SCAN (filename,'.') !IF (dot /= 0) filename = filename(1:dot-1) command = 'gen2shp ' // filename(1:len_trim(filename)) // ' lines < tratti.gen' res = SYSTEM(command) !create dbf command='txt2dbf -C10 -R20.7 -R20.7 -R20.7 -R20.7 -I10 -I10 -R15.7 -R15.2 & -R15.2 -I10 -d ; tratti.csv ' // filename(1:len_trim(filename)) // '.dbf' res = SYSTEM(command) !------------------------------------------------------------------------------ !(1.3) delete intermediate files !------------------------------------------------------------------------------ res = SYSTEM('del tratti.gen') res = SYSTEM('del tratti.csv') RETURN END SUBROUTINE ExportShape END MODULE RiverDrainage